source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/LNP_vaccine/'

lnp_meta <- read_csv('mission/LNP_vaccine/muscle_lnp.meta.csv')

lnp_meta <- lnp_meta |>
  mutate(.cell = ...1,
         prefix = str_remove(.cell, '_[ATCG]{16}.+'),
         sample = case_match(prefix, 'LNP1' ~ 'LNP-2h-1',
                             'LNP2' ~ 'LNP-2h-2',
                             'M-LNP' ~ 'LNP-16h-1',
                             'M_LNP' ~ 'LNP-16h-2',
                             'M_PBS' ~ 'PBS'),
         .keep = 'unused')

mex_lnp <- list.files(proj.nm, 'PBS|LNP', full.names = T, recursive = T) |>
  map(Read10X_h5, .progress = T)

mex_name <- list.files(proj.nm, 'PBS|LNP', recursive = T) |>
  str_remove_all('.+MS_|.h5')

mex_lnp <- mex_lnp |>
  map2(mex_name, add_name_field)

mex_lnp5 <- RowMergeSparseMatrices(mex_lnp[[1]], mex_lnp[-1])

mex_lnp5[1:5,1:5]

sobj <- mex_lnp5 |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj <- sobj |> PercentageFeatureSet('^mt-', col.name = 'mito_ratio')

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- sobj |>
  filter(mito_ratio < 10)

sobj <- sobj |>
  left_join(lnp_meta)

Idents(sobj) <- sobj$sample

sobj <- sobj |>
  filter(!is.na(sample))

sobj |>
  ggplot(aes(group, fill = treat)) + geom_bar()

sobj <- sobj |>
  NormalizeData()

# save RDS ----------
sobj |> write_rds('mission/LNP_vaccine/muscle_LNP.rds')

sobj <- read_rds('mission/LNP_vaccine/muscle_LNP.rds')

# 2h LNP vs PBS --------
lnp2h_frac_logfc <- sobj |>
  filter(str_detect(sample, 'PBS|2h')) |>
  discov_frac_change(group = treat, subtype = annotation,
                     var.1 = LNP, var.2 = PBS)

lnp2h_frac_logfc |>
  ggplot(aes(log2fc_frac, subtype, fill = type)) +
  geom_col() +
  theme_bw() +
  labs(x = 'Proportion change (log2FC)', y = 'Cell type',
       title = 'Cell proportional change in muscle tissue 2hpi of LNP',
       subtitle = 'GSE239574 (Kim et al., Nat Comm 2024)') +
  scale_fill_manual(values = c('royalblue','red2','grey'))

# 16h LNP vs PBS ----------
sobj <- sobj |>
  mutate(group = str_remove(sample, '.\\d$') |> make.names())

lnp16h_frac_logfc <- sobj |>
  filter(group != 'LNP.2h') |>
  discov_frac_change(group = group, subtype = annotation,
                     var.1 = LNP.16h, var.2 = PBS)

lnp16h_frac_logfc |>
  ggplot(aes(log2fc_frac, subtype, fill = type)) +
  geom_col() +
  theme_bw() +
  labs(x = 'Proportion change (log2FC)', y = 'Cell type',
       title = 'Cell proportional change in muscle tissue 16hpi of LNP',
       subtitle = 'GSE239574 (Kim et al., Nat Comm 2024)') +
  scale_fill_manual(values = c('royalblue','red2','grey'))

## cytokine ----------
lnp16_deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'annotation', group.by = 'group',
                       ident.1 = 'LNP.16h', ident.2 = 'PBS')

lnp2_deg <- sobj |>
  FindMarkersAcrossVar(split.by = 'annotation', group.by = 'group',
                       ident.1 = 'LNP.2h', ident.2 = 'PBS')

cytok_gene <- query_uniprot_keyword('KW-0202', species = 'mouse')
chemt_gene <- query_uniprot_keyword('KW-0145', species = 'mouse')
infla_gene <- query_uniprot_keyword('KW-0395', species = 'mouse')

chemokine <- cytok_gene |>
  filter(symbol %in% chemt_gene$symbol)

inflakine <- cytok_gene |>
  filter(symbol %in% infla_gene$symbol)

lnp16_deg |>
  filter(gene %in% chemokine$symbol) |>
  ggplot(aes(gene, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  theme_bw() +
  scale_fill_distiller(palette = 'RdYlBu') +
  rotate_x_text() +
  labs(y = 'Cell type',
       title = 'Chemokine expression in muscle 16 hpi LNP')

lnp16_deg |>
  filter(gene %in% chemokine$symbol) |>
  tidyplot(gene, cluster, color = avg_log2FC, width = 100) |>
  add_heatmap() +
  scale_fill_distiller(palette = 'RdYlBu') +
  labs(y = 'Cell type',
       title = 'Chemokine expression in muscle 16 hpi LNP')

publish_source_plot('chemokine.16hpi.lnp.heatmap', width = 140, height = 70)

lnp16_deg |>
  filter(gene %in% chemokine$symbol, p_val_adj < .05) |>
  tidyplot(gene, cluster, color = avg_log2FC) |>
  add_heatmap() +
  scale_fill_distiller(palette = 'RdYlBu') +
  labs(y = 'Cell type',
       title = 'Chemokine expression in muscle 16 hpi LNP (p<0.05)')

publish_source_plot('chemokine.16hpi.lnp.sig.heatmap', width = 140, height = 70)

lnp16_deg |>
  filter(gene %in% inflakine$symbol) |>
  ggplot(aes(gene, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  theme_bw() +
  scale_fill_distiller(palette = 'RdYlBu') +
  rotate_x_text() +
  labs(y = 'Cell type',
       title = 'Cytokine expression in muscle 16 hpi LNP')

lnp16_deg |>
  filter(gene %in% 'Ifnb1')

lnp16_deg |> write_source_csv('LNP16vsPBS.celltype.deg')

lcf_cytk <- read_csv('mission/LNP_vaccine/lcf.cytokine.csv') |>
  mutate(gene = str_to_title(gene))

lnp16_deg |>
  filter(gene %in% lcf_cytk$gene) |>
  ggplot(aes(gene, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  theme_bw() +
  scale_fill_distiller(palette = 'RdYlBu') +
  rotate_x_text() +
  labs(y = 'Cell type',
       title = 'Cytokine expression in muscle 16 hpi LNP')

lnp2_deg |>
  filter(gene %in% lcf_cytk$gene) |>
  ggplot(aes(gene, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21) +
  theme_bw() +
  scale_fill_distiller(palette = 'RdYlBu') +
  rotate_x_text() +
  labs(y = 'Cell type',
       title = 'Cytokine expression in muscle 2 hpi LNP')

lnp16_deg |>
  filter(gene %in% lcf_cytk$gene, p_val_adj < .05) |>
  tidyplot(gene, cluster, color = avg_log2FC) |>
  add_heatmap() +
  scale_fill_distiller(palette = 'RdYlBu') +
  labs(y = 'Cell type',
       title = 'Cytokine expression in muscle 16 hpi LNP (p<0.05)')

publish_source_plot('cytokine.16hpi.lnp.sig.heatmap', width = 140, height = 70)

lnp2_deg |>
  filter(gene %in% lcf_cytk$gene, p_val_adj < .05) |>
  tidyplot(gene, cluster, color = avg_log2FC) |>
  add_heatmap() +
  scale_fill_distiller(palette = 'RdYlBu') +
  labs(y = 'Cell type',
       title = 'Cytokine expression in muscle 2 hpi LNP (p<0.05)')

publish_source_plot('cytokine.2hpi.lnp.sig.heatmap', width = 140, height = 70)

proj.nm <- 'mission/LNP_vaccine/'
# GO GSEA -----
library(clusterProfiler)

lnp16_deg <-
  read_csv('mission/LNP_vaccine/results/LNP16vsPBS.celltype.deg.csv')

gsego_fibro <- lnp16_deg |>
  filter(cluster == 'Fibroblast', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db', keyType = 'SYMBOL')

gsego_fibro <- gsego_fibro |> simplify()

gsego_fibro |> dotplot()

gsego_fibro@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  distinct(p.adjust, .keep_all = T) |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES', y = 'Pathway',
       title = 'GO BP pathway upregulated in muscle fibroblast 16hpi LNP') +
  theme(plot.title.position = 'plot')

gsego_endot <- lnp16_deg |>
  filter(cluster == 'Endothelial', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(ont = 'ALL', OrgDb = 'org.Mm.eg.db', keyType = 'SYMBOL')

gsego_endot <- gsego_endot |> simplify()

gsego_endot@result |>
  filter(ONTOLOGY == 'BP', NES > 0) |>
  distinct(p.adjust, .keep_all = T) |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES', y = 'Pathway',
       title = 'GO BP pathway upregulated in muscle endothlial cells 16hpi LNP') +
  theme(plot.title.position = 'plot')

# KEGG GSEA --------
lnp16_deg <- lnp16_deg$gene |>
  unique() |>
  bitr(fromType = 'SYMBOL', toType = 'ENTREZID',
       OrgDb = 'org.Mm.eg.db') |>
  as_tibble() |>
  dplyr::rename(gene = SYMBOL) |>
  right_join(lnp16_deg)

gsekegg_fibro <- lnp16_deg |>
  filter(cluster == 'Fibroblast', p_val_adj < .05, !is.na(ENTREZID)) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T) |>
  gseKEGG(organism = 'mmu', eps = 0, use_internal_data = T)

gsekegg_fibro |> dotplot()

gsekegg_fibro@result |>
  filter(NES > 0) |>
  #distinct(p.adjust, .keep_all = T) |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES', y = 'Pathway',
       title = 'KEGG pathway upregulated in muscle fibroblast 16hpi LNP') +
  theme(plot.title.position = 'plot')

gsekegg_endot <- lnp16_deg |>
  filter(cluster == 'Endothelial', p_val_adj < .05, !is.na(ENTREZID)) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T) |>
  gseKEGG(organism = 'mmu', eps = 0)

gsekegg_endot@result |>
  filter(NES > 0) |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES', y = 'Pathway',
       title = 'KEGG pathway upregulated in muscle endothelial cells 16hpi LNP') +
  theme(plot.title.position = 'plot')

gsekegg_monoc <- lnp16_deg |>
  filter(cluster == 'Monocyte', p_val_adj < .05, !is.na(ENTREZID)) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T) |>
  gseKEGG(organism = 'mmu', eps = 0)

gsekegg_monoc@result |>
  filter(NES > 0) |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES', y = 'Pathway',
       title = 'KEGG pathway upregulated in muscle monocytes 16hpi LNP') +
  theme(plot.title.position = 'plot')

### myeloid --------
myelo_16h_deg <- sobj |>
  filter(str_detect(annotation, 'Mono|Neutr|Macro')) |>
  FindMarkers(group.by = 'group',
              ident.1 = 'LNP.16h', ident.2 = 'PBS') |>
  as_tibble(rownames = 'gene')

gsekegg_myelo <- lnp16_deg |>
  distinct(gene, ENTREZID) |>
  right_join(myelo_16h_deg) |>
  filter(p_val_adj < .05, !is.na(ENTREZID)) |>
  pull(avg_log2FC, name = ENTREZID) |>
  sort(decreasing = T) |>
  gseKEGG(organism = 'mmu', eps = 0)

gsekegg_myelo@result |>
  filter(NES > 0) |>
  plot_enrichment(metric = NES) +
  labs(x = 'NES', y = 'Pathway',
       title = 'KEGG pathway upregulated in muscle myeloid cells 16hpi LNP') +
  theme(plot.title.position = 'plot')

# TLR2/4 ----------
sobj |>
  filter(time == '2hr') |>
  DotPlot(c('Tlr2','Tlr4','Myd88'), cols = 'RdYlBu', group.by = 'annotation')

g1 <- last_plot()

g1 + labs(x = 'Gene', y = 'Cell type',
          title = 'Tlr2/4 expression in muscle tissue 2 hpi')

sobj |>
  filter(time == '2hr') |>
  DotPlot(c(tlrs, 'Myd88'), cols = 'RdYlBu', group.by = 'annotation')+
  labs(x = 'Gene', y = 'Cell type',
          title = 'TLRs expression in muscle tissue 2 hpi') +
  RotatedAxis()

tlrs <- sobj |> rownames() |> str_subset('Tlr')
